ergodic.m
Calculates the equilibrium transition rule for (C,N) and its ergodic distribution. If called without Pi, nPrime, nPrimeIndex, minN, maxN, and nCstates in memory, it creates them by calling bellman.m and then creates a plot of the transition rule.
Contents
Call bellman.m to calculate the equilibrium transition rule if required
if ~(exist('Pi','var')&&exist('nPrime','var')&&exist('nPrimeIndex','var')&&exist('minN','var')&&exist('maxN','var')&&exist('nCstates','var')) bellman equilibriumPlot=1; end
Current plot held Current plot held Current plot held Current plot held

Create the transition matrix.
%Stack the columns of nPrime and nPrimeIndex. The relevant Markov chain has as many states as elements of this vector. nprimeIndexVec=nPrimeIndex(:); nprimeVec=nPrime(:); nNstates=maxN-minN+1; %Set up a matrix with indicators for the choice of nprime given (C,N). %This can be big, but it has only ones and zeros. Therefore, use a sparse matrix. nPrimeMask=sparse(ones(nCstates*nNstates,1)*(minN:1:maxN)==nprimeVec*ones(1,nNstates)); %Create the transition matrix for the Markov chain governing (C,N). %This can create quite a large matrix, so a given machine might run out of memory. The try-catch-end construct handles %any resulting errors by swapping as many variables to disk as possible and then attempting the calculation again. try PiCN=nPrimeMask*sparse(kron(eye(nNstates),ones(1,nCstates))); %Create the (deterministic) transition matrix for (C',N)->(C',N'). temp=repmat(Pi,nNstates,nNstates); PiCN=PiCN.*temp; %Multiply each nCstates x nCstates block of PiCN with Pi (elementwise) to create the transition matrix for (C,N) -> (C',N') catch disp('Swapping to disk'); save tempfile save tempfile2 nPrimeMask nNstates nCstates Pi clear load tempfile2 PiCN=nPrimeMask*sparse(kron(eye(nNstates),ones(1,nCstates))); temp=repmat(Pi,nNstates,nNstates); PiCN=PiCN.*temp; clear temp load tempfile end
Calculate the ergodic distribution.
For this, we calculate the eigenvector associated with PiCN's unique unit eigenvalue. Since this is a pretty large matrix and the eigenvector calculation is approximate, some of its elements can be negative. We handle this by setting them to zero and then verifying that the resulting vector's error in the eigenvalue equation is not substantially worse than the original vector's error.
opts.disp=0; %Turn off the display output from eigs. [v,d,flag]=eigs(PiCN',1,'lm',opts); if flag~=0 error('Ergodic distribution calculation failed.'); end v=v/sum(v); %Set any negative elements to zero. if any(v<0) vv=v; vv(vv<0)=0; vv=vv/sum(vv); %Determine whether or not the new value %of v improves on the old one. d1=sum(abs(v'*PiCN-v')); d2=sum(abs(vv'*PiCN-vv')); if d1<d2/10; %The old eigenvector is an order of magnitude better than the new one. error('Calculated ergodic distribution has negative probability values.'); end v=vv; end
Plot the equilibrium transition rule if this is a test run.
This plot displays N' as a function of (C,N). The number's shading indicates its approximate likelihood in the ergodic distribution.
if exist('equilibriumPlot','var') figure(3) set(gcf,'Color',[1 1 1]) CN=[kron(ones(maxN-minN+1,1),omega') kron((minN:1:maxN)',ones(nCstates,1))]; %Set up the plot area. plot(CN(:,1),CN(:,2),'Visible','off') box off; xlim([omega(1) omega(end)]); ylim([minN-0.25 maxN+0.25]); %Assign shadings to different %probabilities. %Find the 'smallest' sets. sortv=sort(v); cdfv=cumsum(sortv); q1index=find(cdfv<=0.005,1,'last'); q2index=find(cdfv<=0.10,1,'last'); q1=sortv(q1index); q2=sortv(q2index); %Pick every tenth state for plotting to avoid clutter. selindx=(1:1:nCstates*(maxN-minN+1)); selindx=find(selindx==10*floor(selindx/10)); CN=CN(selindx,:); nprimeVec=nprimeVec(selindx); vv=v(selindx); %Plot the elements Q2 q2index=find(vv>q2); text(CN(q2index,1),CN(q2index,2),int2str(nprimeVec(q2index)),'FontName','Traditional Arabic','Color',[0 0 0],'FontSize',14,'FontWeight','bold'); %Plot the elements of Q1 not in Q2 q1index=find((vv<=q2).*(vv>q1)); text(CN(q1index,1),CN(q1index,2),int2str(nprimeVec(q1index)),'FontName','Traditional Arabic','Color',[0.25 0.25 0.25],'FontSize',14); %Set the tick marks on the y axis. set(gca,'YTick',minN:1:maxN); ytickstr=int2str(minN); for i=minN+1:maxN ytickstr=[ytickstr '| ' int2str(i)]; end set(gca,'YTickLabel',ytickstr,'FontName','Traditional Arabic','FontSize',16); %Set axis labels xlabel('ln$C_t$','Interpreter','latex','FontSize',22); h=ylabel('$N_t$','Interpreter','latex','FontSize',22); xx=get(h,'Position'); set(h,'Position',[xx(1) xx(2)+0.5]); end
